# -*- coding: utf-8 -*-
###################################################
#    This file is part of blockIMH.
#
#    blockIMH is free software: you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation, either version 3 of the License, or
#    (at your option) any later version.
#
#    blockIMH is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with blockIMH.  If not, see <http://www.gnu.org/licenses/>.
###################################################
#! /usr/bin/env python
source("preliminaries.R")

nbsimu <- 10000
scales <- c(1, 10)

p <- 16
alldf <- list()
melted <- list()
for (indexparameter in 1:3){
    alldf[[indexparameter]] <- data.frame()
    for (scal in scales){
        load(paste("results-probit-p", p, "T", p, "n", nbsimu, "-scale", scal, ".RData", sep = ""))
        df <- data.frame(cbind(basicestimvector[,indexparameter], blockestimvectorRandPermut[,indexparameter],
                               estimRB1vectorRandPermut[,indexparameter], estimRB2vectorRandPermut[,indexparameter]))
        names(df) <- c("basic IMH", "block IMH", "block IMH RB1", "block IMH RB2")
        df$scale <- paste("scale = ", scal, sep = "")
        df$indexsimu <- 1:nbsimu
        alldf[[indexparameter]] <- rbind(alldf[[indexparameter]], df)
    }
    alldf[[indexparameter]]$scale <- factor(alldf[[indexparameter]]$scale, levels = paste("scale = ", scales, sep = ""))
    library(reshape)
    melted[[indexparameter]] <- melt(alldf[[indexparameter]], id = c("indexsimu", "scale"))
}

colourpalette <- hcl(h = seq(120, 0, length = 7) + 150, l = seq(60, 10, length = 7))
finaldataframe <- data.frame()
for (indexparameter in 1:3){
    library(plyr)
    a <- ddply(melted[[indexparameter]], .(scale, variable), colwise(var, c("value")))
    newa <- c()
    for (index in 1:(dim(a)[1])){
        base <- subset(a, variable == "basic IMH" & scale == a[index, 1])
        base <- as.numeric(base[3])
        newa <- c(newa, as.numeric((a[index, 3] - base) / base))
    }
    a$value <- newa
    a <- subset(a, variable != "basic IMH")
    a$variable <- factor(a$variable)
    a$parameter <- paste("parameter", indexparameter)
    finaldataframe <- rbind(finaldataframe, a)
}
finaldataframe$parameter <- factor(finaldataframe$parameter, levels = paste("parameter", 1:3))

g <- ggplot(finaldataframe, aes(x = variable, y = value))
g <- g + geom_bar(aes(size = 1, fill= variable))
g <- g + facet_grid(scale ~ parameter)
g <- g + xlab("") + ylab("") + opts(legend.position = "none")
g <- g + opts(axis.text.x = theme_blank())
g <- g + scale_fill_manual(values = colourpalette)
X11(width = 14, height = 6); print(g)
ggsave(file = "probitComparisonScale.pdf", useDingbats = FALSE)











